%% Fig. 1 Panel A changept algorithm on state level timeseries

%Explanation of variables:
%us_time is a datetime array, which is included in Fig1_policy.mat. 
%At the time of this analysis, data was available from 2/24/2020 to 5/22/2020 (89 days).

%us_e is the Unacast potential encounters metric, which is included in Fig1_policy.mat. 
%county_us_ordered is the equivalence on the county level.


%Output data:
%us_e_change_min is the diagnosed quarantine starts date.
%county_e_change_min and county_e_change_fake_min are its equivalence on
%the county for the real time series and control time series.

%us_e_reopen_change_min is the diagnosed quarantine stops date.
%county_e_reopen_change_min and county_e_reopen_change_fake_min are its equivalence on
%the county for the real time series and control time series.

%us_e_reopen_percent is the date at which potential encounters metric
%reached 30% of its prepandemic value. Its equivalence at the 20%, 40%, and
%50% are included for robustness checks.

us_e_change_min=NaT(51,1);
us_e_reopen_change_min=NaT(51,1);

ind=104; %length of data time series. 
%Though 104 days of data were available when reruning this analysis, the
%initial analysis was already done when there were only 89 days of data
%available.

for i=1:51

    %find 2 changepoints in the time series and make sure they are at least 20 days apart
    tmp=findchangepts(us_e((i-1)*ind+1:(i-1)*ind+89),'MaxNumChanges',2,'MinDistance',20);

    %the first changepoint should correspond to encounter decreasing i.e.
    %going into quarantine
    
    %find the closest local maximum in mobility to the date flagged by the first
    %changepoint. Since the decrease in mobility could be noisy, ensure 
    %maxima are at least 3 days apart to filter out peaks that are too close together.
    [~, i_max,~,p]=findpeaks(us_e((i-1)*ind+tmp(1):...
        (i-1)*ind+40),'MinPeakDistance',3);
    
    %In the continuous period of decreasing mobility, there could be a small
    %blip of local maxima that are clearly still in the period of
    %continuous decrease. The date captured by changept just happens to be to the 
    %left of this blip. To avoid flagging this,look for the most prominent peak, 
    %which would correspond to a clear end of the continuous mobility decrease.
     [~,tmp_peak]=sort(p,'descend');
    
    %if the most prominent peak is the first peak, make sure that it is not a blip
    %for afprementioned reasons. The search for quarantine starts would be before the second
    %peak.
    
    %if the most prominent peak is not the first peak, then it correctly
    %corresponds to when mobility has stabilized. The search for quarantine
    %starts would be between changept and this peak.
     if tmp_peak(1)==1
             [~, i_max]=min(us_e((i-1)*ind+tmp(1):...
         (i-1)*ind+tmp(1)+i_max(tmp_peak(1)+1)-1));
     else
             [~, i_max]=min(us_e((i-1)*ind+tmp(1):...
         (i-1)*ind+tmp(1)+i_max(tmp_peak(1))-1));
     end
      
     %if both changepts happen to both be at the beginning of the time
     %series, then they are both corresponding to the quarantine starts
     %period and the second changept needs to be used. This means
     %findchangepts did not effectively locate a changept for quarantine
     %stops and will be corrected after.
     if length(tmp)==2 && tmp(2)<40
                [~, i_max,~,p]=findpeaks(us_e((i-1)*ind+tmp(2):...
        (i-1)*ind+40));
         [~,tmp_peak]=sort(p,'descend');
        [~, i_max]=min(us_e((i-1)*ind+tmp(1):...
            (i-1)*ind+tmp(2)+i_max(tmp_peak(1))-1));
     
     end
    
     us_e_change_min(i)=us_time(i_max+tmp(1)-1);

    %the second changept should correspond to encounter increasing, i.e. the
    %start of reopening. If the changepts function did not find a second
    %changept or found both in the first half, expand the inputs of the changepts function
    %to 3 to search for this changept. 
     
    if length(tmp)==1 || tmp(end)<40
        tmp=findchangepts(us_e((i-1)*ind+1:(i-1)*ind+89),'MaxNumChanges',3,'MinDistance',10);
    end

    %find the local maximum closest to the changept. No clear evidence to
    %support looking for the most prominent maxima. So the closest local
    %maximum to the changept is used.
    [~, i_max]=findpeaks(us_e((i-1)*ind+tmp(end)-20:...
        (i-1)*ind+tmp(end)));
    [~, i_min]=min(us_e((i-1)*ind+tmp(end)-20+i_max(end)-1:...
        (i-1)*ind+tmp(end)));    
     us_e_reopen_change_min(i)=us_time(i_min-1+tmp(end)-20+i_max(end)-1);

end

%% Fig. 1 Panel A find where encounter rate increase to some percent of precovid
us_e_reopen_percent=NaT(51,1);
percent=0.3;
for i=1:51
    tmp_e=normalize(us_e((i-1)*ind+1:i*ind),'range');
    tmp_base= mean(tmp_e(1:14));%start to 3/8 (inclusive) is pre-COVID
    %average for pre covid baseline
     tmp_t=find(us_time==us_e_change_min(i));
    tmp=find(tmp_e(tmp_t(1)+14:end)>=percent*tmp_base); %54:end
    if ~isempty(tmp)
        us_e_reopen_percent(i)=us_time(tmp(1)+tmp_t(1)+14-1);
    end
end

%% Fig. 1 Panel B changept algorithm on county level timeseries
%essentially the same as the state level algorithm but adapted to work on
%the county level

county_e_change_min=NaT(3054,1);
county_e_reopen_change_min=NaT(3054,1);

county_unacast_fip=unique(table2array(county_us_ordered(:,2))); %extract fip codes for counties

for i=1:length(county_unacast_fip)
    tmp_e=table2array(county_us_ordered(table2array(county_us_ordered(:,2))==county_unacast_fip(i),7));
    tmp_e=tmp_e(1:89);

    tmp=findchangepts(tmp_e,'MaxNumChanges',2,'MinDistance',20);

    if length(tmp)<2
        tmp=findchangepts(tmp_e,'MaxNumChanges',3,'MinDistance',10); 
    end
    
    if length(tmp)>1
        if  tmp(1)<30
            [~, i_max,~,p]=findpeaks(tmp_e(tmp(1):35),'MinPeakDistance',3);   
        else
            [~, i_max,~,p]=findpeaks(tmp_e(tmp(1):tmp(1)+15),'MinPeakDistance',3);
        end

        [~,tmp_peak]=sort(p,'descend');

         if tmp_peak(1)==1 && length(tmp_peak)~=1
                 [~, i_max]=min(tmp_e(tmp(1):tmp(1)+i_max(tmp_peak(1)+1)-1));
         else
                 [~, i_max]=min(tmp_e(tmp(1):tmp(1)+i_max(tmp_peak(1))-1));
         end

         if length(tmp)==2 && tmp(2)<30
                    [~, i_max,~,p]=findpeaks(tmp_e(tmp(2):40));
             [~,tmp_peak]=sort(p,'descend');
            [~, i_max]=min(tmp_e(tmp(1):tmp(2)+i_max(tmp_peak(1))-1));

         end

         county_e_change_min(i)=us_time(i_max+tmp(1)-1);

            %reopen
            [~, i_min]=findpeaks(-1*tmp_e(1:tmp(end)));
            county_e_reopen_change_min(i)=us_time(i_min(end));
     end

 end

%% Fig. 1 Panel B changept algorithm on county level control timeseries 

%essentially same as above but with a different segment of data
%original 89 days 2/24-5/22
%control 89 days 6/1-8/28
county_e_change_fake_min=NaT(3054,1);
county_e_reopen_change_fake_min=NaT(3054,1);

county_unacast_fip=unique(table2array(unacast_county_ordered(:,1)));

for i=1:length(county_unacast_fip)

    tmp_e=table2array(unacast_county_ordered(table2array(unacast_county_ordered(:,1))==county_unacast_fip(i),3));
    tmp_e=tmp_e(99:187);

    tmp=findchangepts(tmp_e,'MaxNumChanges',2,'MinDistance',20);
    if length(tmp)<2
        tmp=findchangepts(tmp_e,'MaxNumChanges',3,'MinDistance',10); 
    end
    if length(tmp)>1
        if  tmp(1)<30 
             [~, i_max,~,p]=findpeaks(tmp_e(tmp(1):40),'MinPeakDistance',3);
        else
            [~, i_max,~,p]=findpeaks(tmp_e(tmp(1):tmp(1)+15),'MinPeakDistance',3);
        end

         [~,tmp_peak]=sort(p,'descend');

         if tmp_peak(1)==1 && length(tmp_peak)~=1
                 [~, i_max]=min(tmp_e(tmp(1):tmp(1)+i_max(tmp_peak(1)+1)-1));
         else
                 [~, i_max]=min(tmp_e(tmp(1):tmp(1)+i_max(tmp_peak(1))-1));
         end

         if length(tmp)==2 && tmp(2)<30
                    [~, i_max,~,p]=findpeaks(tmp_e(tmp(2):40));
             [~,tmp_peak]=sort(p,'descend');
            [~, i_max]=min(tmp_e(tmp(1):tmp(2)+i_max(tmp_peak(1))-1));
         end
         
         county_e_change_fake_min(i)=us_time(i_max+tmp(1)-1+98);
 
        %reopen
        [~, i_min]=findpeaks(-1*tmp_e(1:tmp(end)));
        county_e_reopen_change_fake_min(i)=us_time(i_min(end)+98);
    end
end
